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ABSTRACT 

We describe the development of a flux-limited gray radiation solver for the compressible astrophysics 
code, CASTRO. CASTRO uses an Eulerian grid with block-structured adaptive mesh refinement 
based on a nested hierarchy of logically-rectangular variable-sized grids with simultaneous refinement 
in both space and time. The gray radiation solver is based on a mixed-frame formulation of radiation 
hydrodynamics. In our approach, the system is split into two parts, one part that couples the radiation 
and fluid in a hyperbolic subsystem, and another parabolic part that evolves radiation diffusion and 
source-sink terms. The hyperbolic subsystem is solved explicitly with a high-order Godunov scheme, 
whereas the parabolic part is solved implicitly with a first-order backward Euler method. 
Subject headings: diffusion - hydrodynamics - methods: numerical - radiative transfer 
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1. INTRODUCTION 

In this paper, we present the development of a gray 
radiation solver in our compressible astrophysics code, 
CASTRO. CASTRO uses an Eulerian grid with block- 
structured adaptive mesh refinement (AMR). Our ap- 
proach to AMR is based on a nested hierarchy of 
logically-rectangular variable-sized grids with simultane- 
ous re finement in both spac e and time. In our previous 
paper (jAlmgren et al.l[2010l henceforth Paper I), we de- 
scribe our treatment of hydrodynamics, including gravity 
and nuclear reactions. Here, we describe an algorithm 
for flux-limited gray radiation hydrodynamics based on 
a mixed-frame formulation. 

Many astrophysical phenomena involve radiative pro- 
cesses, which often dominate the energy transport and 
dynamical behavior of the system. Some examples in- 
clude star formation, stellar structure and evolution, ac- 
cretion onto compact objects, and supernovae. Radiation 
hydrodynamics simulations are playing an increasingly 
important role in modeling these astrophysical systems. 

The fundamental equation of radiation t ransfer is a 
six-d imensional integro-differential equation (jPomranind 
Il973f ). which is unfortunately very difficult to solve. Nu- 
merical codes typically solve one or two angular mo- 
ment equations of the transfer equation. One com- 
mon approach is to solve a two-moment system in- 
cludi ng radiation energy density and radiation flux 
(e.g.. iHaves fc NormanI 120031: iHubeny fc Burrow^ 120071: 
iGonzalez et al.ll2007tlSekora fc Stonell2010( ). Thesvstem 
is closed by an approximate expression for the radiation 
pressure. Another popular appro ach is to solve the ra- 
diation ener gy equ ation onl y (e.g.. Turner fc Stonell2001l: 
Haves et al1l2006l: [R rumhol z et al.ll2007l: iSwestv fc Mvral 
20091 : iCommercon et al.ll2011l : lvan der Hoist et al.ll2011f ). 



In this approach, the so-called flux-l imited diffusion 
(FLD) approximation is used for closure (jAlme fc WilsonI 
|1973[) . The two- moment approach is more accurate when 
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the radiation is highly anisotropic and optically-thin, but 
it is computationally more expensive than the FLD ap- 
proach. However, FLD is a very good approximation for 
optically-thick flows. Furthermore, the FLD approach 
can be numerically more robust than the two-moment 
approach for systems in the hyperbolic limit. We have 
adopted the FLD approach in the radiation solver of 
CASTRO. 

CASTRO solves the equations of nonrelativistic ra- 
diation hydrodynamics. Thus, gas quantities, such 
as pressure, temperature and density, are treated as 
frame-independent because the corrections are of order 
0(w^/c^), where v is the gas velocity and c is the speed 
of light. However, radiation quantities, such as radiation 
energy density, radiation flux, and radiation pressure, in 
the comoving frame diff er from those in the labo ratory 
frame by order 0{v/c) (jMihalas fc MihalasI [l999l ) . Ne- 
glecting the 0{v/c) terms potentially leads to erroneous 
results, especially in the dynamic diffusion limit where 
trans port of radiation is dominated by motion of th e 
fluid (ICastod[l972 l : IMihalas fc KleiafTgsl ICastod[200l . 
For numerical codes, some authors chose the comoving 
frame approach in which the radiatio n quantities are 
meas u red in the comoving fram e (e.g.. iTurner fc Stond 
200lt iHaves fc NormaiJ 12001 IGonzalez et all 120071: 



Swestv fc Mvrall2009l : ICommercon et al.ll201l( rwhereas 



others chose the mixed-frame approach in which the ra- 
diation quantities are computed in the laboratory frame 
while the opacities are measured in the comoving frame 
(e.g.. iMihalas fc KleinI [1981 IHubenv fc Burrowsl [20071: 
IKrumholz et al.ll2007t ISekora fc Stonell2010l) . A primary 
weakness of the comoving frame approach is that it does 
not conserve energy, whereas the mixed-frame approach 
is not suitable for systems in which line transport is im- 
portant. In CASTRO, we have chosen the mixed- frame 
approach because it conserves the total energy and is well 
suited for AMR. 

This paper is organized as follows. In §[2]we present the 
governing equations of the mixed-frame gray radiation 
hydrodynamics and the mathematical characteristics of 
the system. In §[3|we describe the single- level integration 
scheme. In §|4|we describe how the integration algorithm 
is extended for AMR. In §[5|we show the scaling behavior 
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of CASTRO with radiation. In § [5] we present results 
from a series of test problems. Finally, we summarize 
the results of the paper in § [71 

2. GRAY RADIATION HYDRODYNAMICS 

2.1. Equations of Gray Radiation Hydrodynamics 

Assuming local thermodynamic equilibrium, the 
mixed-frame frequency-integrated radiation hydrody- 
namics e quations, correct to ( v/c), can be writt en as 
(see e.g., IMihalas fc Kleinlll982l: ILowrie et ai]|1999l ): 



where xr is the Rosseland mean of the sum of the ab- 
sorption and scattering coefficients, and A is the flux lim- 
iter. We adopt the flux l imiter approximation given in 
iLevermore fc Pomraningl ()1981f ) as 
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d{pu) 
dt 



d{pE) 
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| + V.(p«)=0, 



(1) 



The correspond ing radiation pressure tensor is 
(|Levermordll984[ ) 

p(") = i[(l-/)l + (3/-l)nn]i?(«), (11) 



Kp(li)(a7"4 _ ^(o)-|^ ^2'^ where I is the identity tensor of rank 2, n 



V • {pEu+pu) = - cKp{aT'^ - eI°^) 



Xf(-)-F(0), 
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1 dFr 



-Xf{-)-F. 

c 



i?(°)) 
(0) 



(3) 



(4) 



+ np{-){aT^ - E",). (5) 
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Here p, u, p, T, and E are the mass density, velocity, 
pressure, temperature, and total energy per unit mass 
(internal energy, e, plus kinetic energy, u'^/2), respec- 
tively. Er, Fr, and P,. are radiation energy density, ra- 
diation flux, and radiation pressure tensor, respectively. 
Note that here the r subscript denotes radiation. The 
speed of light and radiation constant are denoted by c 
and a, respectively, where a = Aa/c and cr is the Stefan- 
Boltzmann constant. Kp and xf are the Planck mean 
and flux mean interaction coefficients, both in units of 
inverse length. The (0) superscript denotes the comoving 
frame. Radiation quantities {E^, Fr and P^) without the 
(0) superscript are measured in the lab frame. Radiation 
quantities measured in the comoving a nd lab frames are 
related by the Lorentz transformation (jMihalas fc Kleinl 
Il982( ) . It should be noted that absorption and scattering 
coefficients are always computed in the comoving frame 
in the mixed-frame approach. Also note that scattering 
can be included in the flux mean interaction coefficient. 
The whole system is closed by an equation of state for 

the fluid and a relation between and Er'^^ , 

PW =f(°)i?(°), (6) 

where f is the Eddington te nsor in the comoving fr ame. 

In the FLD approximation (|Alme fc Wilsonlll973l ). the 
comoving radiation flux is written in the form of Pick's 
law of diffusion, 

where the diffusion coefficient D is given by 



m 



XR 



(8) 



Vii^r°Vl^-^^°^|, and the Eddington factor / is given by 
f = \ + \^R^. (12) 

We note that in the optically-thick limit both the flux 
limiter A and the Eddington factor / approach 1/3, 
whereas in the optically-thin limit the flux limiter A and 
the Eddington factor / approach and 1, respectively. 

Furthermore we assume that xf = XRj a com- 
mon approximation accurate in optica Uy-thick regions 
' Mihalas fc MihalasI 119990 . Following iKrumholz et al.l 



2003), we keep terms up to 0{v/c), and we drop all 



terms that are insignificant in all following regimes: 
streaming, static diffusion, and dynamic diffusion limits. 
Our radiation hydrodynamics equations now become. 



| + V.(pn) 



0, 



d{pu) 
dt 



d{pE) 
dt 

dEr 



V • [puu) + Vp + WEr = 0, 
V • {pEu + pu) + Am • VEr = 



3-/ 



ErU 



—cKp{aT'^ 
Xu ■ VEr = 



(13) 
(14) 
(15) 

(16) 



cA, 



CKp{aT^ ~ Ef^) -f- V- ( —\IEr 
VXR 

The absorption terms on the right hand side of these 
equations still include the radiation energy density in the 
comoving frame, because this is the the frame in which 
emission and absorption balance as the material becomes 
opaque. The comoving and lab frame quantities are re- 
lated by 



Sl") ^ Er 



+ O(wVc') (17) 
0{v^l^) (18) 
in the framework of the FLD approximation. 



A u 

Er + 2 VEr 

XR c 



2.2. Mathematical Characteristics of the Hyperbolic 
Subsystem 

Radiation hydrodynamics under the assumption of 
FLD is a mixed hyperbolic-parabolic system with stiff 
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source terms. The equations of the hypcrboUc subsys- 
tem are 



ing right eigenvectors are, 



dp 
dt 



djpu) 
dt 



V • (pu) = 
V • (puu) + Vp + X\7Er = 



dt 
dEr 



3-/ 



ErU - Am • \/Er = 0, 



0, 


(19) 


0, 


(20) 


0, 


(21) 


0, 


(22) 



which are obtained by neglecting the terms on the right- 
hand-side of Eqs. [TSHTOl In the hmit of strong equihb- 

rium (i.e., E^^'' « aT^ and xr oo), these right-hand- 
side terms are neghgible and the fuU system becomes 
hyperbohc, governed by Eas.[T5H^ In the more general 
case the hyperbolic subsystem will be solved first as part 
of a time-split discretization. 

In CASTRO, we solve the hyperbolic subsystem with 
a Godunov method, which utilizes a characteristic-based 
Riemann solver. The Godunov method requires that we 
analyze the mathematical characteristics of the hyper- 
bolic subsystem. For simplicity, let us consider the sys- 
tem in one dimension, which can be written in terms of 
primitive variables as. 



where the primitive variables are 



Q = 



and the Jacobian matrix is 



(23) 
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(25) 
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The system is hyperbolic because the Jacobian matrix 
is diagonalizable with four real eigenvalues. In general 
cases, the eigenvalues and eigenvectors are unfortunately 
very complicated. However, when the following relation 
holds, 



3-/ 



A + 1, 



the four eigenvalues are, 

U ~ Cs, U, U, U + Cs, 

where 



Cs = \h- + A + 1) 

VP P 



(26) 



(27) 



(28) 
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and the corresponding left eigenvectors are. 



(0, 
(0, 
(1, 
(0, 



-p/2cs 



0, 
0, 



p/2c. 



l/2c,2, A/2c,2), 

-{X+l)Er/pcJ, Jp/pCs^), 

-l/c,2, -A/c,2), 

l/2c,2, A/2c,2). 



(30) 



These four eigenvectors define the characteristic fields for 
the one-dimensional system. By computing the prod- 
uct of the right eigenvecto rs and the gra dients of their 
corresponding eigenvalues (jLeVeauell2002D . we find that 
the first and fourth fields are genuinely nonlinear corre- 
sponding to either a shock wave or a rarefaction wave. 
The second and third fields are linearly degenerate corre- 
sponding to a contact discontinuity in either gas pressure 
or density. Note that there is also a jump in radiation 
energy density accompanying the jump in gas pressure 
such that the total pressure, ptot = P + XE^, is constant 
across the contact discontinuity. Obviously, in three di- 
mensions, there are two additional linear discontinuities 
for transverse velocities just like the case of pure hydro- 
dynamics. 

It should be noted that Eq. [26] is satisfied in both 
optically-thick and thin limits. Although this condition 
is n ot always satisfied, it is a f airly g ood approximation. 
The iLevermore fc: Pomranin^ ()1981lf fiux limiter that we 
use (Eqs.m&Hni) satisfies 



0.978 < ^—^ - A < 1.05. 



(31) 



Thus, for the purpose of an approximate Riemann solver, 
the approximate eigenvalues and eigenvectors are used. 

2.3. Radiation Diffusion and Source-Sink Terms 

The parabolic part of the system consists of the radia- 
tion diffusion and source-sink terms, which were omit- 
ted from the discussion of the hyperbolic subsystem 
(Eqs.[19]-[22l). 



d{pe) 
dt 



dEr 

~dr 



-c«:p(ar4-£;W) 

— CKp(aT'^ — Er) -\ 
CKp(aT4-S(°)) + 
CKp {aT'^ 



2A— M • VEr 



V • 



XR 
cA, 

XR 



-VEr 



Er) - 2X — U ■ VEr 
XR 



cA, 

XR 



(32) 
(33) 

(34) 

-^Er 

(35)- 



is the specific internal energy. 



The term 



is the radiation modified sound speed. The correspond- 



where e 

CKp(ar^ — Er^') represents the energy exchange in the 
comoving frame between the material and radiation 
through absorption and emission of radiation. The term 
2A(kp/xr)m • V Er is due to the Lorentz transformation 
of radiation energy density. We do not directly solve 
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Eqs. 131] & mi because of our mixed-frame approach. In- 
stead, we solve Eqs. |33] & [35l An implicit treatment 
is usually necessary in order to solve these equations be- 
cause of their stiffness. However, the Lorentz transforma- 
tion term 2A(kp/xr)m • V Er can be treated explicitly in 
many situations because it is of similar order to the term 
\u-\/Er in the hyperbolic subsystem (Eqs.[lT]&[52|), un- 
less the Planck mean is much larger than the Rosseland 
mean. Without scattering, the Planck mean is usually 
larger than the Rosseland mean because the latter gives 
more weight to the lower opacity part of the radiation 
spectrum. However, in many astrophysical phenomena 
(e.g., shock breakout in core-collapse supernovae), elec- 
tron scattering, which contributes to the Rosseland mean 
only, is the dominant source of opacity, and therefore the 
Lorentz transformation term can be neglected in those 



3. SINGLE-LEVEL INTEGRATION ALGORITHM 

For each step at a single level of refinement, the state 
is first evolved using an explicit Godunov method for the 
hyperbolic subsystem (§ 13. Then an implicit update 
for radiation diffusion and source-sink terms is performed 

(§1321). 

It is customary in time-split schemes to denote inter- 
mediate quantities with a fractional time index such as 
n + 1/2, so that, for example, the explicit hyperbolic up- 
date would advance radiation energy density from i?" to 
j^n+i/2 ^j^^ implicit update would then advance it 

from Er^'^^'^ to E'"^^. We are not using this notational 
convention here mainly to avoid confusion in the follow- 
ing section with time-centered quantities constructed at 
the actual intermediate time In § 13.21 where 

we write out the implicit update in detail, we will re- 
fer to the post-hyperbolic intermediate quantities as E~ , 
{pe)~ , etc. 

3.1. Explicit Solver for Hyperbolic subsystem 

The hyperbolic subsystem is treated explicitly. This 
explicit part of our numerical integration algorithm for 
radiation hydrodynamics is very similar to the hydrody- 
namics algorithm presented in Paper I of this series. We 
refer the reader to Paper I for detailed description of the 
integration scheme, which supports a general equation 
of state, self-gravity, and nuclear reactions. Here we will 
only present the parts specific to radiation hydrodynam- 
ics. 

The advection part of the time evolution can be written 
in the form 

. <») 

where U = {p, pu, pE, E^)^ with the superscript T de- 
noting the transpose operation are the conserved vari- 
ables, and F is their flux. The conserved variables are 
defined at cell centers. We predict the primitive vari- 
ables, including p, it, p, pe, Er, from cell centers at time 
i" to edges at time and use an approximate Rie- 

mann solver to construct fluxes, p"+^/^^ on cell faces. 
This algorithm is formally second-order in both space 
and time. The time step is computed using the stan- 
dard CFL condition for explicit methods, with additional 



constraints if additional physics (such as burning) is in- 
cluded. The sound speed used in the computation is now 
the radiation modified sound speed Cg (Eg. \2E\\ . 

3.1.1. Construction of Fluxes 

CASTRO solves the hyperbolic subsystem of radia- 
tion hydrodynamics with an unsplit piecewise parabolic 
method ( PPM) with characterist ic tracing and full corner 
coupling ()Miller fc Colellall2002D . The four major steps 
in the construction of the face-centered fluxes, p'^+^/^^ in 
the case of hydrodynamics have been described in details 
in Paper I. The extension to radiation hydrodynamics 
is straightforward given the characteristic analysis pre- 
sented in ^ 12.21 Thus, we will not repeat the details here. 
We will also omit the procedures for passively advected 
quantities, auxiliary variables, gravity, and reaction, be- 
cause they do not change. 

First, we compute the primitive variables defined as 
Q = {p,u,p, pe, Er,ptot, pe + Er)^ . Note that several 
variables are redundant for various reasons (i.e., effi- 
ciency and safety), which will be explained later. 

Second, we reconstruct parabolic profiles of the prim- 
itive variables within each cell. The total pressure, ptot, 
rather than the gas pressure, p, is used in computing a 
flattening coefficient ()Miller fc Colellall2002D . 

In the third step, we perform characteristic extrapola- 
tions of the primitive variables and obtain the edge val- 
ues of Q at using the eigenvectors of the system 
(Eq. [251 & [50]) and the parabolic proflles of the primitive 
variables. Flattening is applied in this procedure. 

Finally, the fluxes are computed for the edge values 
obtained in the last step using an approximate Rie- 
mann solver, wh ich i s based on th e Riem ann solver of 
IBeh et all (ITOSl) and lColella eTall (119971) . The compu- 
tational procedure is essentially the same as that in Pa- 
per I except that the gas internal energy density, pe, and 
gas pressure, p are now replaced by the total internal 
energy density, pe + Er, and total pressure, p + XEr- 
The Riemann solver computes the Godunov state at 
the interface, which is then used to compute the fluxes, 
{pu, puu + p\, pEu + pu, ((3 — f)/2)Eru)'^ . Recall that 
there are redundant variables in the primitive variables, 
Q, for efficiency and safety. With the total internal en- 
ergy density, pe + Er, a call to the equation of state can 
be avoided. Negative radiation energy density and gas 
pressure can also be avoided in the Godunov state. The 
term. Am • WEr, in Eqs.HH&j^His computed as follows. 



(Am • VEr)tj.k = XijM 



(37) 



,UxA-l/2,j,k + 'U-xA+l/2J,k ■.,E'r,i+i/2,j,k ~ ^r,i-l/2 j",fe > 
^ 2 A^^ ' 

,'^y.i,j~l/2,k + '"■y.ij"+l/2,fc w-E'r.ij"+l/2.fe ~ -E'r,i.j-l/2.fc ^ 

+^ 2 ' 

/^z,ij\fe-l/2 + , , £'r,ij,fc+l/2 ~ Er^ij^k-1/2 x 

^ 2 A'z '. 

where the variables with half-integer index are the Go- 
dunov states. The term, WEr, in Eq. [201 is computed 
in a similar way. 

Depending upon a switch set by the user, the Lorentz 
transformation term, 2A(Kp/xr{)M • VEr, may be in- 
cluded in the explicit update. 
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3.2. Implicit Solver for Radiation Diffusion and 
Source-Sink Terms 

The implicit solver evolves the radiation and g 0jS 0,C- 
cording to Eqs.[33]&|3S] The algorithm uses a first-order 
backward Euler discretization. We note that the Lorentz 
term, 2A(kp / xrI^- ' V-E^, may or may not be included in 
the implicit solver. The advantage of including this term 
here is that it effectively balances emission against ab- 
sorption in the comoving frame as the material becomes 
optically-thick and this coupling becomes stiff. Implicit 
treatment also helps avoid a time step restriction when 
Kp/xR is significantly greater than 1. The disadvantage 
is that it makes the resulting linear systems nonsymmet- 
ric, but this not a major concern in practice. 

The radiation update al gorithm is based on that of 
iHowell fc Greenoughl ()2003[ ). The update from the post- 
hyperbolic state to time t"^^ for Eqs. [33l fc [35] has the 
form 



At 

- E- 

At 

-,n+l 



.n+1 



(38) 



71 + 1 



n+l 



.n+1 



P - (39) 



where = 2A"+1k;^+Vxr+^ and = 

^^Ti+i^^n-i-i^ Here, we use the — superscript to denote 
the state following the explicit update, and the n-\-\ su- 
perscript for the state at Since the velocity does not 
change in the implicit radiation update, we have dropped 
the — superscript for u~ . We solve Eqs. [38l fc l39l itera- 
tively via Newton's method. We define 



pe"+i - pe 



At{-cK^+^ [a{T'+^f - 

+ q^+^u-VE'^+^}, (40) 
Fr = E'^+^ - E- - At {ci<^+^ [a(r"+i)4 - E'^+^] 

~q''+^u ■ \/E';+'^ + V • {d"+^VE;'+^ )} . (41) 

Here, we have dropped the — superscript for p~ because 
the implicit update does not change the mass density. 
The desired solution is for F,, and Fr to both be zero, 
and the Newton update to approach this state is the 
solution to the linear system 



(dFe/dT)^''^ 


{dFjdErY^^ ' 




" ST^k+l) ■ 






(dFr/dT)^''') 


(dFr/dEr-Y''^ 




. SEi''+'^ 







(42) 

Here (5T(^+i) = T'^+hik+i) _ j^n+Uk) ^nd = 
j^n+i,{k+i) _ j^n+i,(k) ^ where the (fc) superscript denotes 
the stage of the Newton iteration. To reduce clutter 
we drop the n + 1 superscript without loss of clarity. 

To solve this system we eliminate the dependency on 
5T by forming the Schur complement, leaving a modified 
diffusion equation for the radiation update 5Ej.. For sim- 
plicity we drop the temperature derivatives of d and q, 
keeping only the temperature dependence of the emission 
and absorption terms. This does not affect the converged 
solution and in practice does not appear to significantly 



degrade the convergence rate. 

After some mathematical manipulation we obtain the 
following diffusion equation, which must be solved for 
each Newton iteration: 



(1 - ry)cKp + 



1 

At 



p{k+i) 



V • [dVEl!'+^^ 



(43) 



= (1 - 77)cKpa (rC^) 
where 

PCv 



+ (1 - vi)qu 
1 

'At 



VE, 



(k+i) 



e: 



rj{pe} 



pe 



PCv 



. d r 
cAt— [kp 



aT'^ - Er)] 



dFr 



OF, 
dT 



(44) 

and Ci, is the specific heat capacity of the matter. 
The iterations are stopped when the maximum of 

\5Ei^^'^'' / Er^'^''^''^'^\ on the computational domain falls 
below a preset tolerance (e.g., 10^®). Note that Eq. l43l 
is rewritten to be in terms of ^^'+^'(''"+1) rather than 

5Er'''^^\ This is done for computational efficiency. After 
one or two Newton iterations, the solution at the previ- 
ous iteration, Er'^^'^'^\ is getting very close to the final 



n+l,(fc) 



as the 



converged solution E"""*"^. Since we use E., 
starting point for the call to the iterative linear solver 
these calls get cheaper for each additional Newton iter- 
ation. If we solved for Se!/''^^'^ instead, we would have 
to change the linear solver tolerance to avoid an unnec- 
essarily accurate and expensive solve for what may be a 
very small correction. 

Each time we solve the diffusion equation for a new it- 
erate TffQ update the gas internal energy den- 
sity as follows, 



pe 
-At{l 



{k+i} 



v) 



ripe'-"' + (1 - 
a(T('=))4 



CKp 



f])pe 



(45) 

qu ■ v£;('=+i)" . 



The new temperature then derives from 

gn-i-i,(fc-l-i) g^jjjj p yjg^ pg^^j ^Yie equation of state. Other 
quantities may or may not be updated: the coefficients 
f^Pi XRj d, and q are also temperature-dependent and 
could have been written with a (fc) superscript. The lim- 

iter A can be recomputed based on the new Er ' 
It is a tradeoff between efficiency and accuracy whether 
to recompute some or all of these at every iteration. Our 
default choice is to recompute all of these at each it- 
eration for better accuracy. However, updating coeffi- 
cients can make the implicit update iteration less stable 
if the coefficients arc not smooth functions of their in- 
puts, in addition to the extra computational costs. But 
since each stage of the Newton iteration (Eq. 23] com- 
bined with Eq. |45)) is itself a conservative solution, the 
implicit algorithm will conserve total energy to the tol- 
erance of the linear solver (which should not be confused 
with the tolerance of the Newton iterations) regardless of 
the number of iterations taken or which coefficients arc 
updated. 

In CASTRO, the linear sy stem Eg. l43l is solved 
using the hypre library (|Falgout fc Yand 120021 : 
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\hvvre Code Pro iec^ 120 111 ). We have developed drivers 
that work with systems in the canonical form 



ae: 



n+l 



E 



6> 

dx'- 



E 



d 
dx^ 



dE] 



n+l 



dx^ 



rhs, 



(46) 



where A are cell-centered coefficients and Bi, Ci, and 
Di are centered at cell faces. For Eq. 231 the Ci coeffi- 
cients are not used. There is some subtlety in the ap- 
propriate averaging o f coeffi cients from cells to faces; see 
iHowell fc Greenoughl ()2003D for further discussion. The 
same canonical form works for Cartesian, cylindrical, and 
spherical coordinates so long as appropriate metric fac- 
tors are included in the coefficients and the rhs. 

With hypre we have a choice between two parallel 
multigrid solvers: SchafFer multigrid (SMG) and PFMG. 
SMC is more robust for difficult problems with strongly- 
varying coefficients, but PFMG is typically more efficient 
and scalable. These solvers work for systems at a uniform 
grid resolution (that is, systems associated with a single 
level of adaptive mesh refinement). For systems coupling 
together more than one refinement level we could use the 
hypre al gebraic multigrid ( AMG) solver, or an FAC-type 
scheme (IMcCormickl 119891) using structured solvers on 



the separate levels. In earlier versions of the AMR algo- 
rithm we required multilevel solvers for conservative cou- 
pling between refinement levels. We have now developed 
a scheme, though, that eliminates the need for a multi- 
level linear solver while still conserving total energy. We 
discuss this in detail in the following sections. Note also 
that use of the Di coefficients deriving from the Lorentz 
term make the diffusion equation nonsymmetric. The 
multigrid solvers mentioned above are designed for use 
with symmetric systems, but good convergence behavior 
can still be obtained by using these solvers as precondi- 
tioners for a Krylov method such as GMRES. 

4. AMR 

CASTRO uses a nested hierarchy of logically- 
rectangular, variable-sized grids with simultaneous re- 
finement in both space and time, as illustrated in Fig- 
ures [T] and [2j One major design objective of the AMR 
algorithm is to preserve the conservation properties of the 
uniform-grid discretization. AMR for hyperbolic equa- 
tions in CASTRO was described in detail in Paper I. 

The explicit update step for radiation hydrodynam- 
ics (§ I3.ip follows the same pattern as other hyperbolic 
equations and so does not increase the complexity of the 
AMR algorithm. We note that the hyperbolic subsys- 
tem (Eqs. [T9l- [22|) is only partially in conservation law 
form. It will not conserve total momentum because it 
does not include an equation analogous to Eq. [5] track- 
ing the radiation momentum. It will conserve total en- 
ergy, though, so long as the divergence terms are differ- 
enced in a conservative manner. These divergence terms 
therefore require AMR reflux operations, as described in 
iBerger fc Colellal (|I989|. The term Xu ■ VEr appears in 
both the gas energy and radiation energy equations with 
opposite signs, so any consistent discretization of these 
terms will conserve total energy. 

The AMR version of the implicit radiation diffusion up- 




FlG. 1. — A properly nested hierarchy of grids. Each grid consists 
of a number of cells. Thick lines represent level boundaries. The 
union of fine grids at level I is contained within the union of coarser 
grids at level I — 1. 



H (7)—^ Level 2 



-<2> 



-H Level 1 



H Level 



Time 

Fig. 2. — One coarse time step for an adaptive run with one base 
level and two refinement levels. The numbers mark the order of 
the steps. 



date is based on IHowell fc Greenoughl (|2003D . but with 
several important differences: The present algorithm is 
fuUy-implieit, not time-centered. The optional multilevel 
linear solve at the beginning of each coarse time step is 
no longer included — this feature was introduced to im- 
prove accuracy but we now consider it unnecessary in 
most cases. Finally, the multilevel linear solve for flux 
synchronization between coarse and fine levels is replaced 
by a new algorithm we call the "deferred sync." These 
changes entirely eliminate the need to compute linear 
system solutions coupling different levels of the AMR 
hierarchy, while not compromising conservation of total 
energy. Performance is significantly improved because 
multilevel linear solvers tend to be more complex and 
expensive than those for single-level systems. 

4.1. AMR Time Step Outline with Deferred Sync 

The AMR time step is defined recursively in terms of 
operations on a level i and its interactions with coarser 
and finer levels. We consider advancing level i from time 
index n to n -I- 1, corresponding to time values t°^^'^ and 
^now,^^ respectively. (Even though levels other than £ 
have executed different numbers of time steps, we will use 
the n-\-l superscript to refer to values at time f™'"'*' on all 
levels involved in the calculation.) The region covered by 
level I is denoted A^, its border is 9A^, and the border of 
the next finer level, projected onto level £, is P(9A^"''^). 

The notation to describe all of this is unavoidably com- 
plex due to the quantities at different times, levels, and 
stages of the update process. In the following outline, we 
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specify the update for the level i and its synchronization 
with finer levels. We include the hyperbolic update and 
the refluxing step associated with it in order to show the 
proper sequence of operations and to contrast the ex- 
plicit reflux with the implicit deferred sync. As in § 13.11 
the hyperbolic conserved state vector is denoted by U , 
but we denote the hyperbolic flux by i^H to distinguish 
it from the radiation flux. For the radiation flux we are 
concerned here only with the diffusion term in the im- 
plicit update, and we denote the associated flux by 
to distinguish it from the complete radiation flux F^. in- 
troduced in § 12.11 

Note that while the flux divergence is needed every- 
where, the fluxes Fh and Fr themselves are stored only 
on the borders between levels. Our code has data struc- 
tures called flux registers designed for this purpose. The 
notation (•) indicates an average of level £ + 1 data in 
space over the fine cells (or cell faces) making up each 
corresponding coarse cell (or face), while ((•)) denotes 
an average of level £ + 1 data in both space and time 
over the coarse (level £) time step. Thus, at the only 
point in the algorithm where this notation is used, 



1 



/ + i(„+l)- 

E 



y+l,m+l\ 



(47) 



where r^'^^ is the number of level £+1 time steps making 
up the level £ time step from n to n-l- 1, and p'^+^'"^+^ is 
the level £+1 flux during the level £+1 time step 171 + I. 

The expression involving VRoflux' that appears as a 
deferred source term in the diffusion equation and ex- 
plicitly for the hyperbolic reflux is called the reflux di- 
vergence, because it takes the form of the divergence of 
a flux difference 5F stored in the flux registers. These 
terms are evaluated only in the coarse cells bordering an 
interface with a flner level. They do not affect fine cells 
at the interface because flux calculations for those cells 
are already the more accurate ones; the VroAux' terms 
represent the corrected fluxes from these fine cells being 
imposed onto the coarse grid. 

Another way to understand the VroAux' terms is to 
consider 5F^ and 6F^ to be energy that has been "mis- 
placed" at the coarse-fine interfaces during the level time 
step, due to the differing flux calculations on the differ- 
ent levels. If the solution were not corrected, this energy 
would be lost and the system would not be conservative. 
Instead, the VroHux' terms re-introduce the missing en- 
ergy into the system. For the explicit hyperbolic flux this 
is done explicitly; for the radiation flux it contributes to 
the right hand side for the implicit update, so as not to 
impose a new stability constraint on the size of the time 
step. 

The hyperbolic state vector U includes the radiation 
and fluid energies updated in the implicit update step, 
but there is no ambiguity because the operations and 
their associated fluxes are completely distinct. The re- 
fluxing update to U is written as an update, with U 
appearing on both sides of the equation, because other- 
wise we would need additional notation to indicate the 
pre- and post-reflux states. Averaging down from fine to 
coarse levels is also written as an update. The meaning 
of the rest of the pseudocode below should be reasonably 
clear in context: 



If (£ < €inax) and (regrid requested from base level £) 
then 

For f g{C-l,...,£} do 

• Determine new grid layout for level £' + 1. 

• Interpolate data to new grids from level £'. 

• Copy data on intersection with old level £' + 1. 
Enddo 

Endif 

Level Time Step, level £: 

Explicit Hyperbolic Update for J7^ "+i 

(see Paper I) 

• Set ^^^"+1 for hyperbolic fluxes 

on (P(9A^+i), £ < £,„ax) and on ((9A^ i > 0) 
Implicit Diffusion Update 

e:}+^ - E- 



V • (rf"+iv£;;'+i) 

(At'VAr+l) VRefiux • (-^-F^^''") 

on A^ 



71+1 — 

End Implicit Diffusion Update 



on 



A^ 



• F 



I,n+1 
R 



on (P(aA^+i), £ < £„iax) and on {dK\ £ > 0) 
• Advance levels ^ H- 1 , . . . , ^max recursively. 



on P(9A^+i), £<£„ 
+1 

on P(9A^+i), £<£„ 



End Level Time Step 

If {£ < 4iax) then 

(synchronization/refluxing between levels £ and £ + 1) 

. U'-+' U'--+' - (Ar+i) VroHux • (5F^H^^'"+^) 

on A^-P(A^+i) 

on P(A^+i) 



Endif 



,n+l 



(U 



e+l,n+l\ 



(end synchronization / refluxing) . 



The advantages of the deferred sync algorithm are that 
it eliminates the need for a separate linear solve for syn- 
chronization and eliminates the need for solving multi- 
level linear systems entirely. It does, however, add com- 
plications of its own to the AMR implementation. One 
is that there is no longer any point within the time step 
cycle when the field variables are fully synchronized. At 
the end of a coarse time step all levels have reached the 
same point in time, but if we want to actually compute 
the energy budget and confirm conservation we have to 
include the contributions from deferred fluxes stored in 
flux registers. These will not be re-introduced into the 
state until the next time step. The end of a coarse time 
step is also the natural time for checkpoint /restart oper- 
ations, so to reproduce the saved state of the system we 
now have to include the deferred fluxes in checkpoints as 
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Fig. 3. — Weak scaling behavior of CASTRO on Hopper at 
NERSC. Average wall eloek times per coarse time step are shown 
for simulations with 1 (circle), 2 (X symbol), and 3 total grid levels 
(triangle). The numbers of cells that are advanced in one coarse 
time step increase by a factor of three and seven, for the two- and 
three-level simulations, respectively. 

well. 

Regridding becomes an issue as well. The adaptive 
algorithm periodically re-evaluates the refinement crite- 
ria for each level and may change the layout of refined 
grids. For grids at level £ the criteria are evaluated at 
level i — 1, and this happens between level £ — 1 time 
steps. The interpolation operations between coarse and 
fine field variables are conservative. With the deferred 
sync, though, there will also be fluxes stored around the 
edges of level £ at the time that level may be changed. 
There is no straightforward way to transform these stored 
fluxes so that they coincide with the new mesh layout. 
Instead, what we have is an old set of flux registers that 
may now overlap with level £ as well as with level £ — 1, 
and if the grid layout changes enough may even overlap 
with other levels both finer and coarser. The deferred 
sync idea still applies, but the implementation becomes 
more complicated than the pseudocode above suggests. 
Portions of the stored flux at the interface between levels 
£ and £ — 1 may be reintroduced into the level £ advance 
or the level £ + 1 advance, and so on, not just into level 
£-1. 

5. PARALLEL PERFORMANCE 

CASTRO is implemented within the BoxLib frame- 
work for parallel structured-grid AMR applications (Pa- 
per I and references therein). In BoxLib, paralleliza- 
tion is based upon either hybrid OpenMP-MPI or pure 
MPI. Because the hypre library does not fully support 
OpenMP yet, we use the pure MPI approach for the ra- 
diation hydrodynamics solver. For more information on 
software design and parallelization, we refer the reader 
to Paper I. Here we show the scaling behavior of the 
radiation hydrodynamics solver in CASTRO. 

A weak scaling study has been carried out on HopperQ, 
a petascale Cray XE6 supercomputer at the National 
Energy Research Scientific Computing Center. We have 
performed a series of three-dimensional simulations with 
1, 2 and 3 total levels on various numbers of cores. For 

* The Hopper supercomputer was named after Grace Hopper, a 
pioneer in computer science and the developer of the first compiler. 



the convenience of comparison, each run has one grid of 
64^ cells at each level on each core. A refinement factor of 
2 is used in the multi-level simulations. Thus, the level 
1 and 2 grids occupy 12.5% and 1.5625% of the whole 
volume, respectively. A point explosion like the one in 
? l6.9l is replicated on each core. The fine grids are placed 
at the center of the local domain of each core. Figure [3] 
shows that CASTRO has very good scaling behavior up 
to 32768 cores. For the single-level simulations, the aver- 
age wall clock time per coarse time step increases by only 
17% from 8 cores to 32768 cores. Because of subcycling 
in time for simulations with multiple levels, one coarse 
time step consists of one step on the coarse level and 
two steps on the next fine level, and the two fine steps 
might also consist of even finer steps, recursively. Thus, 
the number of cells that are advanced in one coarse time 
step increases by a factor of three or seven, for the two- 
and three- level simulations, respectively (Fig. [5]). Our 
results also show that the overhead introduced by AMR 
is modest. For example, on 4096 cores, the average wall 
clock times per coarse time step for the two- and three- 
level simulations are 3.6 and 8.6 times more than that 
of the single-level simulation. This corresponds to an 
overhead of 19% and 22%, respectively. On 32768 cores, 
the average wall clock times per coarse time step for the 
two- and three-level simulations are 3.8 and 9.5 times 
more than that of the single-level simulation. This cor- 
responds to an overhead of 25% and 36%, respectively. It 
should be noted that single-level simulations with equiv- 
alent uniform grids would cost ~ 4.5 and 30 times more 
time than the corresponding two- and three-level simula- 
tions we have run even if the single-level simulations are 
assumed to scale perfectly. We also note that in this se- 
ries of simulations about half of the time is spent on the 
implicit evolution of the parabolic part of the system. 

6. TEST PROBLEMS 

In this section we present detailed tests of the code 
demonstrating its ability to handle a wide range of ra- 
diation hydrodynamics problems. Note that not every 
term is included in every test so that the algorithms for 
these terms can be tested separately. We test the ra- 
diation source-sink term in isolation in the approach to 
thermal equilibrium test (§ 16. 1|) . For a nonequilibrium 
Marshak wave problem (§ 16. 2[) . the simulation involves 
the parabolic subsystem only. In § 16.31 we assess the 
accuracy of the AMR algorithms for radiation diffusion 
using a thermal wave test, and perform a convergence 
study. The system in a radiation front test (§ 16. 4p is in 
the optically-thin streaming limit, whereas the system in 
a shock tube problem (§ 16. 5|) is in the limit of strong 
equilibrium with almost no diffusion. The full radiation 
hydrodynamics system is included in a nonequilibrium 
radiative shock problem (§ 16. 6p and the advection of a 
radiation pulse problem (i; l6.7p . In 16.81 we demonstrate 
the ability of CASTRO to maintain a static equilibrium 
of the gas and radiation pressures. In a radiative blast 
wave test 16. 9p . wc compare the results of simulations 
in ID spherical, 2D cylindrical (r and z), and 3D Carte- 
sian coordinates. Finally, we demonstrate the ability of 
CASTRO to handle a large Lorentz transformation term 
in another radiative blast wave test (? I6.10p . 

A CFL number of 0.8 is used for these tests unless 
stated otherwise or a fixed time step is used. The refine- 
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Fig. 4. — Evolution of internal energy density of gas for calcu- 
lations of the approach to thermal equilibrium. Numerical results 
are shown in symbols, whereas the analytic solutions are shown 
in solid lines. Two cases with an initial internal energy density 
of lO^^ergcm"^ (upper solid line and squares) and 10 erg cm~^ 
(lower solid line and circles) are shown. The time step is fixed at 
At = 10"" s. 

ment factor is 2 for all AMR runs. The relative tolerance 
for the Newton iterations in the implicit update is 10~^ 
for all runs. The Lorentz transformation term is handled 
explicitly except in the second radiative blast wave test 
(§«. 

6.1. Approach to Thermal Equilibrium 

This test introduced bv lTurner &: Stond ()2001[ ) has of- 
ten been used to test the ability of a code to handle 
the source-sink term, cKp(aT^ — Er). The test con- 
sists of a static uniform field of gas and radiation. The 
gas has a density of p = 10~^gcm~'^, a Planck mean 
absorption coefficient of Kp = 4 x 10~^cm~^, a mean 
molecular weight of /x = 0.6, and an adiabatic index 
of 7 = 5/3. The initial radiation energy density is 
Er = lO^^ergcm^'^ corresponding to a temperature of 
~ 3.39 X 10^ K. Two cases with an initial internal energy 
density of 10^ erg cm^"^ and 10^°ergcm~^, respectively, 
have been studied. The gas temperatures are ^ 4.81 K 
and ~ 4.81 x 10^ K for the two cases, respectively. The 
time step is chosen as At = 10~^-'^s. The evolution of the 
system will bring the gas and radiation into a thermal 
equilibrium. The radiation energy density will hardly 
change because the energy exchange between the gas and 
radiation is only a small fraction of the radiation energy. 
Therefore an analytic solution can be calculated by solv- 
ing the following ordinary differential equation 



djpe) 
dt 



'CKp{aT'^ — Er), 



(48) 



where Er is assumed to be constant. The results of this 
test of the approach to thermal equilibrium are shown 
in Figure |4l The agreement with the analytic solution 
is good, especially in the first case. For the second case, 
the cooling in the numerical calculation is initially slower 
than that in the analytic soluti on for the first few step s 
in agreement with the results of lTurner fc Stoiii ()2001[ ). 
The slow cooling is a result of the backward Euler method 
and a relatively large time step at the beginning of the 
decay of the gas temperature. 



Fig. 5. — Nonequilibrium Marshak wave. Numerical results are 
shown in lines, whereas the analytic solutions are shown in circle 
symbols. We show the dimensionless radiation energy density u at 
T = 0.01 {solid line) and t = 0.3 (dash-dot line), and the dimen- 
sionless gas energy density v at r = 0.01 (dashed line) and t = 0.3 
(dotted line). Here t is the dimensionless time. 

6.2. Nonequilibrium Marshak Wave 

In this test, we simulate the nonequilibrium Mar- 
shak wave problem in one dimension. Initially half of 
the space, z > 0, consists of a static uniform zero- 
temperature gas and no radiation. A constant radiation 
flux i^inc is incident on the surface at z = 0. The gas is 
not allowed to move in this idealized test. Thus the gas 
and radiation are coupled only through the source-sink 
term an d the systeni is goy erned by Eqs. [33] & [351 with 
u ~ 0. iSu fc OlsonI (|1996[ ) obtained analytic solutions 
for the problem under special assumptions. The matter 
is assumed to be gray with Kp = xr, and its volumet- 
ric heat capacity at constant volume is assumed to be 
Cy = aT^, where T is the gas temperature and a is a 
parameter. 

We have run this test with 4a/a = 0.1 and no flux 
limiter (i.e., A = 1/3). Figure [5] shows the numerical 
results of the dimensionless ra diation energy den sity and 
gas energy density defined bv iPomranind (|l979l ) 



u{x, t) 
v(^x, t) 



( 4acK\ 

Er{z,t) 
^ inc 

fc\ ( aT^{z,t) 



\ a 

id 



4/ V ^ir 



(49) 
(50) 

(51) 

(52) 



The computational domain of < x < 5-\/3 is covered 
by 128 uniform cells. The dimensionless time step is 
chosen to be At = 3 x 10~^. The numerical results are in 

ood agreement with the analytic results of iSu fc OlsonI 

19961) . 

6.3. Thermal Wave 
In this test, we simulate a thermal wave 



(jZel'Dovich fc Raizerlll967l) in three dimensions and use 
this test to assess the accuracy of the AMR algorithms 
for radiation diffusion. Suppose that a large amount of 
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Fig. 6. — Temperature and error at t = 0.006 s for the deferred sync run of the thermal wave test. The error is computed as the difference 
between the simulation results and the analytic solution. We show a 2D slice at x = 0. Also shown is the grid structure of the adaptive 
mesh. The physical domain in each dimension is (—200, 200) cm for each panel. 



energy is deposited into a small volume as the internal 
energy of matter. The matter then starts to radiate 
and transfer most of its energy to radiation. We assume 
that the Planck mean absorption coefficient is large 
enough so that the matter and the radiation are in 
thermal equilibrium. The heat is transported out of the 
initial hot spot because of the nonlinear radiation heat 
conduction. As a result, a thermal wave develops. As 
the matter cools down, the matter gains most of the 
energy again. Initially the thermal wave speed is much 
higher than the sound speed. Assuming that there is 
no fluid motion and the matter contains most of the 
energy, there is an analyt ic solution for this problem 

(iZel'Dovich fc Raized[T967l) 

This test is adapted from lHowell fc Greenoughl ()2003[ ). 
The computational region in this test is three- 
dimensional with a domain of (—200, 200) cm in each 
dimension. Initially, the spherical hot spot has an 
energy of 3 x 10^ erg within a radius of 3.125 cm, 
whereas the ambient medium has a very low tempera- 
ture of 10~^K. The volumetric heat capacity is pc^ = 
0.05 ergK~^ cm""^. The Planck mean absorption coeffi- 
cient is Kp — lO^cm"^, whereas the Rosseland mean is 
XR = 10"3(y/iK)V2cm-i. In this test, the hydrody- 
namics is turned off, and there is no flux limiter (i.e., 
A = 1/3). We have performed three simulations, a non- 
AMR run with 128^ uniform cells, an AMR run using 
the deferred sync algorithm (§ 14. I J. and an AMR run 
using the multilevel algorithm of iHowell fc Greenoughl 
()2003D . The two AMR runs use 2 total levels with an 
effective resolution of 128'^ cells on the finer level. A cell 
is tagged for refinement if its temperature satisfies both 
VT > OAT/ Ax and T > 10"^ K, where Ax is the size 
of the cell. For the non-AMR run, we use a time step 
At = 1.03"~^ X 5 X 10~^^s for step n, whereas for the 
two AMR runs, we use At = 1.0609"-^ x 1.015 x lO^^^s 
for step n on the coarser level. We have chosen the time 
steps for the following reasons. First, as it expands, the 



thermal wave is rapidly decelerated. Thus, a fixed time 
step would not be optimal. Instead, we let the time 
step grow over time. Second, the numerical error de- 
pends on both time step At and cell size Ax. In this 
test, we want to assess the accuracy of the AMR runs in 
comparison with a non-AMR run. Therefore, we make 
the time step on the finer level of the AMR runs to be 
roughly the same as that of the non-AMR run. Figure [S] 
shows a 2D slice at a; = of the deferred sync run at 
t = 0.006 s for temperature and the difference between 
the numerical and analytic results. It is shown that the 
largest errors {'^ 5-10 K) occur near the thermal wave 
front where the slope of the temperature profile is ex- 
tremely steep. At the center, the absolute error is only 
0.6 K, whereas the relative error is 0.7%. It takes 451 
coarse time steps for the AMR runs to reach t = 0.006 s. 
The finer grids in the AMR runs occupy 0.20, 2.3, and 
32% of the volume at steps 1, 226, and 451, respectively. 
Thus, the benefit of AMR is obvious. To quantitatively 
measure the accuracy of the results, we compute the nor- 
malized ii-norm error for temperature, j ^. \T-n.ij,k — 
T^{tn, x,,yj,Zk)\AVtj^k/ Eij.fc Ti,{tn, Xi,yj,Zk)AV,j,k , 
where Tn^ij^k and Ta(tn, Xi, yj, Zk) are the numerical and 
analytic results for cell {i,j,k), respectively, and AVi,j^k 
is the volume of cell (i, j, k). At t = 0.006 s, the Li-norm 
errors are 0.00792, 0.00833, and 0.00847, for the non- 
AMR, deferred sync, multilevel sync runs, respectively. 
The results show that our AMR algorithms have only 
increased the error by 5%-7%. 

We have also performed a series of AMR simulations 
using the deferred sync to check the convergence behav- 
ior of the code. Besides the AMR run using the deferred 
sync that has been presented (Fig. two lower resolu- 
tion runs and one higher resolution run have been per- 
formed (Tabled]). In afl four runs, there are two total 
AMR levels. In the convergence study, when the cell size 
changes by a factor of 2 from one run to another, we 
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TABLE 1 

Ll-NORM ERRORS AND CONVERGENCE RATE 
FOR THE THERMAL WAVE TEST PROBLEM. 

Four 3D AMR runs with various 

RESOLUTIONS ARE SHOWN. THERE ARE TWO 
AMR LEVELS IN EACH OF THESE RUNS. 



Resolution^ 


Li Error 


Convergence Rate 


16, 32 


0.0706 




32, 64 


0.0245 


1.5 


64, 128 


0.00833 


1.6 


128, 256 


0.00268 


1.6 



^ Number of cells across the width of the do- 
main at each of the two AMR levels 



change the time step by a factor of 4. Table [T] shows the 
Li-norm errors and convergence rate aX t = 0.006 s. In 
this study, the order of accuracy with respect to Ax is 
'--^ 1.6. Because our imphcit scheme is first-order in time 
and second-order in the spatial discretization of the dif- 
fusion term, the expected convergence rate for a smooth 
flow is second-order with respect to Aa; when At/Aa;^ is 
fixed. However, the temperature profile of the thermal 
wave problem has a very steep slope near the front and 
its second derivative is discontinuous there. Hence, it 
is not surprising that the achieved order of accuracy is 
lower than 2. 

6.4. Optically-Thin Streaming of Radiation Front 

In this problem, we test our code in the optically-thin 
streaming limit. This is the opposite limit from diffu- 
sion. The computational domain of this test is a one- 
dimensional region of < x < 100 cm, covered by 128 
uniform cells. Initially, the region a; < 10 cm is filled 
with radiation, Ej. = lergcm"'^, whereas the region of 
X > 10 cm has Er ~ 10~^° erg cm~^. The left and right 
boundaries are held at Er = 1 and 10~^° erg cm"'', re- 
spectively, during the calculation. The hyperbolic up- 
date (§ 13. ip is switched off in this test. The Planck 
mean absorption coefficient is set to zero. We have 
studied two cases with Rosseland means of xr = 10"^ 
and 10~^cm~^, respectively. Thus the optical depths 
of the whole domain are 10"^ and 10~^, respectively. 
For each case, we have performed two calculations, one 
with the time step set to At = Ax /{2c) and one with 
At = A.t/(20c). Physically, the radiation front is ex- 
pected to propagate at the speed of light, but the diffu- 
sion approximation does not naturally support a propa- 
gating front at any speed. Only the flux limiter permits 
the code to approximate the correct solution. In the 
numerical results the radiation front moves at approxi- 
mately the speed of light with spreading due to diffusion 
(Figure [7]). Better results are obtained for smaller time 
steps. 

6.5. Shock Tube Problem In Strong Equilibrium Regime 

In this test, the one-dimensional numerical region (0 < 
X < 100 cm) initially consists of two constant states: 
= 10-5gcm-3, Tl = 1.5 X lO^K, ul = and 
= 10-5gcm-3, Tfl = 3 X lO^K, ur = 0, where L 
stands for the left state, and R the right state. The ini- 
tial discontinuity is at x = 50 cm. The gas is assumed 
to be ideal with an adiabatic index of 7 — 5/3 and a 
mean molecular weight of /i = 1. Initially, the radia- 
tion is assumed to be in thermal equilibrium with the 
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Fig. 7. — Optically- Thin Streaming of Radiation Front. We 
show radiation energy density at t = 2 X 10~^ s for two cases: 
(^) Xr = 10~*cm~ (upper panel); (b) xr = 10~^cm~^ (lower 
panel). Results with a time step of At = Aa;/(2c) are shown in 
solid lines, whereas those with a time step of At = Aa;/(20c) in 
dotted lines. The vertical lines at a; = 70 cm indicate the expected 
position of the radiation front. 
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Fic;. 8. — Shock tube at i = 10^^ s. Numerical results from a full 
radiation hydrodynamics calculation with 128 cells arc shown in 
symbols. Results from a hydrodynamics calculation with 128 cells 
and an EOS for ideal gas plus radiation are shown in sold lines for 
comparison. We show mass density (p), velocity (?;), total pressure 
(ptot), and radiation energy density (-Er). The total pressure is the 
sum of gas pressure and radiation pressure. All quantities are in 
cgs units. 

gas (i.e., Er = aT'^). The interaction coefficients are 
set to Kp = lO^cm"^ and xr = 10® cm"^. Thus, due 
to the huge opacities, the system is close to the limit 
of strong equilibrium with no diffusion (i.e., Er « aT* 
and Xr oo), and is essentially governed by Eqs. [19]- 
[221 The parameters of this test problem are chosen such 
that neither radiation pressure nor gas pressure can be 
ignored. The numerical results are shown in Figure [5) 
Also shown are the results from a pure hydrodynamics 
calculation with an equation of state (EOS) for ideal gas 
plus radiation. As we expected, the results from the ra- 
diation hydrodynamics calculation are almost identical 
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to those of the pure hydrodynamics calculation. We also 
note that our scheme is stable without using small time 
steps even th ough the system is in the "dynamic dif- 
fusion " limit (jMihalas fc MihalasllTOOl iKrumholz et alJ 
I2007D with Tu/c ^ 10^, where r is the optical depth of 
the system. 

6.6. Nonequilihrium Radiative Shock 

Radiation can modify the structure of a shock be- 
cause it diffuses and because it interacts with matter. 
Analyti c estimates of radiative sho c k structures can be 
found in IZerPovich fc Raized (I1967D : iMihalas fc MihalasI 
(|1999[ ). Recentlv lLowrie fc Edwardsl (Eooli " have found a 
semi-analytic exact solution of radiative shocks with gray 
nonequilihrium diffusion. In this section, we present our 
numerical results for tw o shock strengths and co mpare 
them to the solutions of iLowrie fc Ed wards' (2008 ^). The 
first problem is the Mach 2 shock in iLowric fc Edward j 
dlOOS^). This is a subcritical shock in which the pre-shock 
matter is preheated by diffused radiation to a tempera- 
ture that is lower than the temperature in the relaxation 
region far downstream. There is an embedded hydrody- 
namic shock that causes a jump in density and raises the 
temperature of the matter above the final downstream 
temperature. The matter temperature behind the em- 
bedded hydrodynamic shock then cools down in the re- 
laxation region. The seco nd problem is the Mach 5 super- 
critical shock in lLowrie~fc Edwards (2008). In the case of 
a supercritical shock, the matter temperature just before 
the hydrodynamic shock equals the downstream temper- 
ature. Most of the precursor region of a supercritical 
shock is close to equilibrium. Note that in both cases 
there is no jump in radiation energy density, otherwise 
the radiation flux would be infinite, which is unphysical. 

The Mach 2 shock problem is simulated in a one- 
dimensional region of —1000 cm < .t < 500 cm consisting 
of two constant states: pL = 5.45887 x 10^^^ g cm~^, 



and PR 
= 1.02987 X 



Tl = lOOK, UL = 2.35435 x lO^cms-i 
1.24794 X 10-12 gcm-^ Tr = 207.757K, ur 
lO^cms^i, where L stands for the left state, and R the 
right state. The initial discontinuity is at x = 0. The 
gas is assumed to be ideal with an adiabatic index of 
7 = 5/3 and a mean molecular weight ol p = 1 . Initially, 
the radiation is assumed to be in thermal equilibrium 
with the gas (i.e., = aT^). The Planck and Rosse- 
land coefficients are set to Kp = 3.92664 x 10~^ cm"^ and 
Xr = 0.848902 cm~i, respectively. The boundaries are 
held at fixed values. Two refinement levels (three total 
levels) are used with the finest resolution at Ax « 2.9 cm. 
In this test, a cell is tagged for refinement if the normal- 
ized seco nd derivative of eit her density or temperature 
given by ()Frvxell et al.|[2000l ) 



\ui+2 - Ui\ + \ui - Ui-2\ + e(|Ui+2| + 2|ui| -f- \Ui^2\) 

(53) 

is greater than 0.8. Here, e is a parameter set to 0.02 in 
this test, and u denotes either density or temperature. 
The initial conditions are set according to the pre-shock 
and post-shock states of the M2 shock. After a brief 
period of adjustment, the system settles down to a steady 
shock. The simulation is stopped at t = 0.05 s. The 
results are shown in Figure |9l The agreement with the 
analytic solution is excellent, except that the numerical 




X (cm) 

Fig. 9. — Density and temperature profiles for Mach 2 subcritical 
shock. Numerical results are shown in symbols. We show density 
{circle), gas temperature {plus sign), and radiation temperature 
{square). Here, radiation temperature is defined as {Er /a)^^'^ . The 
analytic solution is shown in solid lines. We show only part of 
the region near the hydro-shock. The numerical results have been 
shifted by +10 cm in space to compensate for the discrepancy in 
shock position caused by the initial setup. 




X (cm) 

Fig. 10. — Density and temperature profiles for Mach 5 su- 
percritical shock. Numerical results are shown in symbols. We 
show density (circle), gas temperature [plus sign), and radiation 
temperature {square). Here radiation temperature is defined as 
{Er/a)^^^. Also shown are the analytic results for density {red 
solid line), gas temperature {blue solid line), and radiation tem- 
perature {green solid line). We show only part of the region near 
the hydro-shock. The inset shows a blow-up of the spike in tem- 
perature. The numerical results have been shifted by —205 cm in 
space to compensate for the discrepancy in shock position caused 
by the initial setup. 
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results in the figure had to be shifted by 10 cm in space to 
match the analytic shock position. This discrepancy in 
shock position is due to the initial transient phase as the 
initial state relaxes to the correct steady-state profile. 
No flux limiter is u sed in this c alculation because the 
analytic solution of ILowrie fc Edwards (.2008:) assumes 
A= 1/3. 

The setup of the Mach 5 shock problem is similar to 
that of the Mach 2 shock problem. The computational 
domain in this test is —4000 cm < a; < 2000 cm. Four 
refinement levels (five total levels) are used with the 
finest resolution at Ax ~ 1.5 cm. The refinement cri- 
teria are the same as in the Mach 2 shock test (Eq. [55)1 . 
The initial left and right states are pL — 5.45887 x 
10-13 gcni-3, Tl = 100 K, UL = 5.88588 x lO^cms^^ 
and PR = 1.96405 x lO'^^ g ^.^^^-3^ j,^ ^ 855.720K, 
ur = 1.63592 X lO^cms--^, respectively. Again, the sys- 
tem settles down to a steady shock after a brief period 
of adjustment. Figure [TU] shows the results at i = 0.04 s. 
The results including the narrow spike in temperature 
are in good agreement with the analytic solution. 

6.7. Advecting Radiation Pulse 

In this test, introduced bv lKrumholz et all (|2007f ). we 
simulate the advection of a radiation pulse by the motion 
of gas. The initial temperature and density profiles are 




Fig. 11. — Density, velocity and temperature profiles at t = 
4.8 X 10~^ s for the test of advecting radiation pulse. The results 
of three runs are shown. The difference is so small that it is nearly 
invisible to the eye. 



r = To + (ri-To)exp --^ 



2w' 



Tn 



P - 



afi 



Tl 3i? V T 



(54) 
(55) 



where To = 10^ K, Ti = 2 x 10^ K, po = 1.2 gcm-^, 
w — 24 cm, the mean molecular weight of the gas is 
p, = 2.33, and R is the ideal gas constant. Initially the 
radiation is assumed to be in thermal equilibrium with 
the gas. We also assume that the radiation pressure is 
Er/3 (i.e., A = 1/3). If there were no radiation diffusion, 
the system would be in an equilibrium with the gas pres- 
sure and radiation pressure balancing each other. The 
interaction coefficients are set proportional to density as 
'^p = Xr = lOOcm^g^i X p. Because of radiation diffu- 
sion, the balance is lost and the gas moves. The purpose 
of this test is to assess the ability of the code to handle 
radiation being advected by gas. We will solve the prob- 
lem numerically in different laboratory frames. Then we 
can compare the case in which the gas is initially at rest 
to the case in which the gas initially moves at a constant 
velocity. 

We have performed three runs for compariso n with 
each o ther and against the results of fKrumhol z et al.l 
(j2007| ). The computational domain in all three runs 
is a one-dimensional region of —512 cm < a; < 512 cm 
with periodic boundaries. The simulations are stopped 
at t = 4.8 X 10"^ s. The velocity in the first run is 
initially zero everywhere, whereas in the second run it 
is lO^cms"^ everywhere. The numerical grid in these 
two runs consists of 512 uniform cells. The third run is 
a high-resolution calculation of the first case with 4096 
uniform cells. Figure [Tl] shows the density, velocity, and 
temperature profiles for all three runs. The results of 
the advected run have been shifted in space for com- 
parison. The profiles from these three runs are almost 
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Fig. 12. — Relative errors in density {solid lines) and gas tem- 
perature [dashed lines) in the test of advecting radiation pulse. 
The difference between the low-resolution unadvected and high- 
resolution unadvected runs is shown in panel (a), whereas the dif- 
ference between the advected and unadvected runs in panel (b). 



identical demonstrating the accuracy of our scheme in 
this test. We also show the relative difference in Fig- 
ure [T21 The relative error of the low-resolution unad- 
vected run with respect to the high-resolution run is 
computed as (high - low) /high, and the relative error of 
the low-resolution advected run with respect to the low- 
resolution unadvected run is computed as (unadvected 
- advected) / unadvected. Note that in Figure [T^] the 
results of the high-resolution run have been restricted 
to the low-resolution grid by averaging for comparison. 
The figure shows that the difference between the low- 
resolution and high-resolution runs is less than 0.05% 
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Fig. 13. — Structure of the magnitude of velocity at i = 10~*s 
in the static equiUbrium test. The unit of velocity is cm s~^. 

everywhere, and the difference between the advected and 
unadvected runs is less than 0.0016% everywhere. The 
re lative error in ou r runs is 1000 times smaller than that 
of iKrumholz et all (j2007| ). Wc also note that the asym- 
metry in the difference between the advected and unad- 
vected cases is expected because our upwinding scheme 
breaks the symmetry slightly in the advected run. The 
linear solver can also break the symmetry slightly. 

6.8. Static Equilibrium 

In the test of advecting radiation pulse f? l6.7p . if there 
were no radiation diffusion, the system would be in a 
static equilibrium with the gas and radiation pressures 
balancing each other. In this test, we set both interac- 
tion coefficients to a very high value of 10^" cm^ g~^ x p. 
Thus almost no diffusion can happen. We have per- 
formed a calculation on a two-dimensional Cartesian grid 
of —512 cm < X < 512 cm and —512 cm < y < 512 cm 
with 512 uniform cells in each direction. The initial ve- 
locity is zero everywhere. For the initial setup, the coor- 
dinate X in Eq.[54lis replaced by r = y^c^ + y^. Figure [T51 
shows the velocity profile at t = lO^'* s. The maximal ve- 
locity at that time is '--^ 4 x 10~^cms~"'^. Note that the 
sound speed in this problem is between 2.5 x 10^ and 
6.1 x 10^ cm s^^. Such a small gas velocity indicates that 
CASTRO can maintain a perfect static equilibrium in 
multiple dimensions because of the way radiation pres- 
sure and gas pressure are coupled in the Riemann solver. 

6.9. Radiative Blast Wave: Case 1 

In this test, a large amount of energy is deposited into a 
small region. This results in a spherical blast wave, which 
is somewhat similar to the Scdov- Taylor blast wave in 
hydrodynamics, except here we include radiation. There 
is no analytic solution for this problem. 

We use this problem to test the implementation of the 
geometric factors in CASTRO. We have performed simu- 
lations in ID spherical, 2D cylindrical (r and z), and 3D 
Cartesian coordinates. Initially, the gas is at rest with a 
density of 5 x 10~^gcm~'^. The initial temperature for 
both the gas and radiation is set to 10'^ K, except for a re- 
gion inside a sphere with a radius of 2 x 10^^ cm centered 
at the origin where temperature is set to 10^ K. The gas 
is assumed to be ideal with an adiabatic index of 7 = 5/3 
and a mean molecular weight of /i = 1. The Planck and 
Rosseland coefficients are set to Kp = 2 x 10~^^ cm^^ and 
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Fig. 14. — Density, radial velocity and temperature profiles at 
t = 10^ s for the first radiative blast wave problem. Here, radial 
velocity is the velocity in spherical radial direction. We show the 
results of the ID spherical run {red), 2D cylindrical run (green), 
3D Cartesian run (blue), and the high-resolution ID spherical run 
(black). The bottom panel shows both the gas temperature (solid 
lines) and radiation temperature (dashed lines). The inset shows 
a blow-up of the spike in gas temperature, where the ID, 2D, 3D, 
and the high-resolution runs arc shown in red circles, green plus 
signs, blue triangles, and black lines. 

Xr = 2 X 10~^°cm~^. These simulations were run with 
two refinement levels (three total levels) with a cell size of 
^ 3.9 X 10^^ cm at the finest level. In these simulations, 
at each but the finest level a numerical cell is tagged for 
refinement if it satisfies either p > 5.01 x 10~^ g cm~'^ or 
T > 9 X 10^ K. Furthermore, the intermediate level is 
slightly larger than the finest level due to proper nest- 
ing (see Paper I). The computational domain for the ID 
spherical run is < r < lO^'* cm. The computational 
domain for the 2D cylindrical run is < r < 10^"* cm 
and —10^^ cm < z < 10^^ cm. For the 3D run, the com- 
putational domain is (—10^^, 10^''') cm in each direction. 
A CFL number of 0.6 is used for these simulations, and 
the initial time step is shrunk by a factor of 100 to al- 
low the point explosion to develop. Note that we have 
chosen these initial conditions so that different regimes 
can be explored by the test. The blast wave starts with 
almost all the energy in radiation. At the end of the sim- 
ulations, approximately one third of the energy is in the 
gas, and the gas pressure just behind the shock is about 
twice the radiation pressure. Furthermore, the gas and 
radiation are not in equilibrium near the shock due to 
the low Planck mean opacity. 

Since there is no analytic solution for the problem, we 
have also run a high-resolution ID simulation with a cell 
size of ^ 4.9 x 10^° cm for comparison. Figure [Ml shows 
the radial profiles of density, radial velocity, gas temper- 
ature, and radiation temperature ai t ~ 10^ s for these 
runs. The radial profiles of the 2D and 3D results are 
computed by mapping each cell into its corresponding 
radial bin and averaging. The width of the bins is cho- 
sen to be the cell size at the finest refinement level. The 
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Fig. 15. — Snapshot of the 2D cylindrical simulation at t = 10 s for the first radiative blast wave test. Profiles of density, radial velocity, 
gas temperature, and radiation temperature are shown in panels (a), (b), (c), and (d), respectively. Here radial velocity is the velocity in 
spherical radial direction, not the cylindrical radial direction. The physical domain in each dimension is (—10^*, 10^*) cm for each panel. 
Also shown is the grid structure of the adaptive mesh. Levels 0, 1, and 2 are shown in black, white, and purple, respectively. 




Fig. 16. — 2D slice at 2 = of the 3D Cartesian simulation at t = 10 s for the first radiative blast wave test. Profiles of density, radial 
velocity, gas temperature, and radiation temperature are shown in panels (a), (b), (c), and (d), respectively. Here radial velocity is the 
velocity in spherical radial direction. The physical domain in each dimension is (—10^*, 10^^) cm for each panel. Also shown is the grid 
structure of the adaptive mesh. Levels 0, 1, and 2 are shown in black, white, and purple, respectively. 
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Fig. 17. — Density, radial velocity and temperature profiles at 
t = 10^ s for the second radiative blast wave problem. The results 
of the 2D cylindrical simulation are shown in symbols {black plus 
signs). We do not show the radiation temperature because it is 
almost identical to the gas temperature. Here, radial velocity is 
the velocity in spherical radial direction, not the cylindrical radial 
direction. Also shown are the results of the high-resolution ID 
spherical simulation (red solid lines). 

results from runs in three different coordinates are in ex- 
ceUent agreement with each other, and they agree with 
those of the high-resolution simulation. A snapshot of 
the structure of density, radial velocity, gas temperature, 
and radiation temperature at t = 10^ s for the 2D simu- 
lation is shown in Figure [15] A 2D slice at z = of the 
3D simulation is shown in Figure 1161 The multi-D re- 
sults show good agreement with each other, and they are 
spherically symmetric except for some minor asymmetry 
in velocity at the low density region near the center. This 
asymmetry is, in part, due to the representation of the 
initial hot sphere in non-spherical coordinates. Another 
source of the asymmetry in the 2D run is the coordi- 
nate singularity at the longitudinal axis (r = 0) of the 
cylindrical coordinates. It should also be noted that the 
AMR gridding algorithm in CASTRO does not necessar- 
ily preserve symmetry even if the layout of the numerical 
cells that are tagged for refinement is symmetric. Last, 
but perhaps not least, the linear solver does not preserve 
perfect symmetry either. 

6.10. Radiative Blast Wave: Case 2 

This test is similar to the first radiative blast wave test 
(§ 16. 9|) except that the Planck mean coefficient is set to 
Kp = 2 X 10^ cm~^ . Thus the ratio of the two coefficients 
is Kp/xR = 1000. Because of such a large ratio of the 
two mean opacities, the Lorentz transformation term has 
to be treated implicitly in this test, otherwise the time 
step would have to be much smaller for stability reasons. 

We have performed a 2D cylindrical simulation with 
two refinement levels (three total levels) and a cell size 
of ^ 3.9 X 10^^ cm at the finest level on a computational 
domain of < r < 10^^ cm and —10^* cm < z < lO'^* cm. 
In this 2D AMR run, a numerical cell is tagged for re- 



finement if it satisfies either p > 5.01 x 10~^gcm~^ or 
r > 9 X 10^ K. We have also performed a high-resolution 
ID simulation with a cell size of ^ 4.9 x 10^° cm for com- 
parison. A CFL number of 0.6 is used for the simulations, 
and the initial time step is shrunk by a factor of 100 to 
allow the point explosion to develop. Figure [17] shows the 
radial profiles of density, radial velocity, gas temperature 
and radiation temperature at t = 10® s for the two runs. 
The results show that CASTRO can handle a very large 
Lorentz transformation term implicitly without having 
to decrease the time step. 

7. SUMMARY 

We have developed a new radiation hydrodynamics 
solver in our compressible astrophysics code, CASTRO. 
We use the mixed-frame approach and adopt the FLD 
assumption. The solver uses a second-order explicit Go- 
dunov method for the hyperbolic part of the system and 
a first-order backward Euler method for the parabolic 
part. We have also presented the mathematical char- 
acteristics of the hyperbolic subsystem. The eigenvalues 
and eigenvectors of the system we have obtained are used 
to construct the Riemann solver in our Godunov scheme, 
and could also be useful for other characteristic-based 
schemes. 

We have demonstrated the capability of CASTRO 
to address a wide range of radiation hydrodynamics 
problems by extensive testing. There are a number 
of other radiation hydrodynam ics codes. Some recen t 
examples i nclude ZEUS-2D (iTurner fc Stond 120011) 
ZEU S-MP ([Haves et al .l I2006D . Orion (iKrumholz et al l 
[200l. H ERACLES (iGonzalez et al.^ I2007D V2p 
([Swestv fc M vra 2001. Nike (fSekora fc Stond [2010h . 
RAMSES ([Commercon et all l201l[) . and CRASH 
([van der Hoist et al.ll2011[ ). Here we compare CASTRO 
with these other codes. 

• A major advantage of CASTRO is its efficiency due 
to the use of AMR and combined with good scal- 
ing behavior on up to 32768 cores. Among the 
other radiation hydrodynamics codes listed above, 
only Orion, RAMSES, and CRASH are three- 
dimensional AMR codes. ZEUS-MP and HERA- 
CLES are three-dimensional codes without AMR. 
ZEUS-2D and V2D are two-dimensional codes. 
Nike is a one-dimensional code. 

• A unique strength of our code is that the hyper- 
bolic solver in CASTRO is an unsplit version of the 
PPM method that avoids spurious noise caused by 
dimensional splitting (see Paper I for an example 
of the advantage of unsplit methods over dimen- 
sional splitting methods). For the hyperbolic part, 
Orion, HERACLES, Nike, RAMSES, and CRASH 
use Godunov methods, whereas ZEUS-2D, ZEUS- 
MP, and V2D use the ZEUS type of algorithms. 
The latter are faster but somewhat less accurate 
than high-order Godunov methods. However, for a 
radiation hydrodynamics code, the cost of the hy- 
perbolic solver is no longer a concern, because the 
implicit parabolic solver is more expensive. For 
multigroup radiation hydrodynamics, the implicit 
parabolic solver would be even more expensive and 
the hyperbolic solver would be essentially free. 
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• Our scheme is very robust even for dynamic dif- 
fusion. By solving the entire hyperbolic subsys- 
tem in one Riemann solver, the scheme can avoid 
the operator splitting errors that appear in many 
other numerical schemes (sec e.g.. iKrumholz et alJ 
[200I . We note that RAMSES and CRASH also 
couple both radiation and matter in their Riemann 
solvers. 

• CASTRO is based on a mixed-frame formulation, 
as is Orion and Nike. A main advantage of the 
mixed-frame approach is that it conserves the total 
energy, whereas its main disadvantage is that it 
is of limited use for line transport. However, line 
transport cannot be treated by a gray radiation 
solver regardless of its choice of frame. 

• CASTRO, ZEUS-2D, ZEUS-MP, Orion, V2D, 
RAMSES, and CRASH have adopted the FLD ap- 
proach, whereas HERACLES and Nike are based 
on the two-moment approach. Thus CASTRO car- 
ries the limitations of FLD, such as poor accuracy 
for optically-thin flows. However, FLD is compu- 
tational cheaper than the two-moment approach. 
Moreover, for multigroup radiation, FLD is much 
less memory intensive than the two-moment ap- 
proach. 

The current implementation uses a gray approximation 
based on the frequency-integrated formulation of the ra- 
diation hydrodynamics equations. We note that numer- 
ical codes exist for multigroup flux-limited radiation hy- 



drodynamics (e.g.. iBurrows et al.l [20071 : iSwestv fc Mvral 
I2009t lvan der Hoist et al.ll201lD . It is straightforward to 
extend our scheme to multigroup radiation because the 
radiation pressure in the hyperbolic subsystem is still a 
frequency-integrated quantity for multigroup radiation. 
A multigroup neutrino-radiation hydrodynamics solver 
is also currently under development. 

Further details on CASTRO can be found iii the CAS- 
TRO User Guide (IC ASTRO User Guide! [201I . 
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